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\ ABSTRACT 

■ We present a study of the pre-protostellar core L1544. A series of self-consistent, 
' three-dimensional continuum radiative transfer models are constructed. The outputs 

t — , of these models are convolved with appropriate telescope beam responses, and com- 

^sO • pared with existing SCUBA data. The resulting comparison allows us to constrain the 

,— 1 | structure of L1544. We find that the source is well-fit by a prolate spheroid, having 

an ellipsoidal power-law density distribution of index m ~ 2 (1.75 < m < 2.25) in 

■ to at least r ~ 1600 AU. For r < 1600 AU, the data are consistent with either an 

extension of the power law to smaller radii, or a flattened (Bonner-Ebert like) density 

• distribtion. Furthermore, we find an optical depth along the short axis at 1300/im 

P-f of Ti3oo :S hort = 5 x 10~ 3 (2 x 10~ 3 < Ti 30 o, short < 8 x 10~ 3 ), a central luminosity 

q ! L» = (< 10~ 3 L Q ), a long axis diameter D = 0.1 pc (0.08 < D(pc) < 0.16 ; 

16000 < D(AU) < 32000), an axis ratio q = 2 (1.7 < q < 2.5), and an external ISRF 
! defined by Mathis, Mezger, & Panagia (1983) to within 50 per cent. The outer diam- 

C$ • eter and axis ratio may each be somewhat larger due to potential on-source chopping 

in the observations, and the projection of the long axis onto the plane of the sky. 
. i-H ! While these results are similar to those inferred directly from observations or sphcr- 

■ ical modeling due to the source transparency at submillimeter wavelengths, we infer 
5-h ' a smaller size, lower mass, and higher optical depth / column density, exposed to a 

stronger external radiation field than previously assumed. Finally, we find that both 
the spectral energy distribution (SED) and surface brightness distribution are neces- 
sary to constrain the source properties in this way, and even a modest variation in \ 2 
can significantly alter the fit quality. 

Key words: stars: formation - infrared: stars - ISM: clouds. 



1 INTRODUCTION 

The structure of bound, starless cores (e.g., Myers, Linke, 
& Benson 1983; Benson & Myers 1989) is not well known. 
This structure is important, as it should drive the ensuing 
evolution of the core, and give insight into the initial condi- 
tions of star formation free from the disruption of a central 
protostar. Ward-Thompson, Scott, Hills, & Andre (1994) 
detected millimeter and submillimeter continuum radiation 
from dust in these regions, and dubbed them preprotostellar 
cores (PPCs). 

Following this detection, significant work has contin- 
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ued toward determining the conditions associated with these 
regions. Ward-Thompson, Scott, Hills, & Andre (1994) in- 
ferred that the density distributions do not follow isother- 
mal sphere distributions (n(r) oc r~ 2 ) over a range of cores, 
a result supported by various other studies (e.g. Andre, 
Ward- Thompson, & Motte 1996; Ward-Thompson, Motte, 
& Andre 1999; Bacmann et al. 2000). These studies utilized 
semi-analytic techniques to infer the dust density distribu- 
tion. Similar techniques such as colour temperature varia- 
tion, were used by Ward- Thompson, Andre, & Kirk (2002) 
(along with ISO flux data) to infer the temperature distribu- 
tion and to place constraints on the luminosity of the central 
source. While many of these semi-analytic approaches yield 
physically meaningful information, the specifics of the ap- 
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proach can affect the reliability of the conclusions (Doty & 
Palotti 2002). 

Submillimeter (submm) data for many PPCs were ob- 
tained from the Submillimeter Common-User Bolometric 
Array (SCUBA) on the James Clerk Maxwell Telescope 
(JCMT) by Shirley et al. 2000. In an analysis of these data, 
Evans et al. (2001) assumed spherical symmetry and uti- 
lized a self-consistent dust radiative transfer model to infer 
the properties and structure of their source sample. They 
found that Bonner-Ebert (BE) spheres provided good fits 
to the data, with LI 544 requiring the highest central den- 
sity, n c = 10 6 cm~ 3 , in fact. However, a power-law density 
distribution could not be ruled out when a self-consistent 
dust temperature variation was included. This result, while 
intriguing, relies upon the assumed underlying source ge- 
ometry. Unfortunately, the observed isophotes for L1544 are 
decidedly non-spherical, leading to potential uncertainty in 
these conclusions. 

In this paper we report on a study comparing existing 
data with detailed, three-dimensional, continuum radiative 
transfer models to constrain the conditions associated with 
L1544. In section two, we describe L1544 and previously in- 
ferred properties. We briefly describe the model and range of 
parameter space considered in section three. In section four, 
we compare the models with observations in constraining 
the source properties. Finally, we conclude in section five. 



2 L1544 

L1544 is located in Taurus at a distance of 140 pc (Elias 
1978). The source has been well- observed. There is no 
reported IRAS source associated with the core (Ward- 
Thompson, Scott, Hills, & Andre 1994). However, submil- 
limeter data have been obtained at the JCMT at 450, 850, 
and 1100 jj,m with a single element by Ward-Thompson, 
Scott, Hills, & Andre (1994), and using SCUBA (Shirley 
ct al. 2000). Further submillimeter data were taken at 1.3- 
mm using IRAM (Ward- Thompson, Scott, Hills, & Andre 
1994; Ward-Thompson, Motte, & Andre 1999). L1544 has 
also been studied with the Infrared Space Observatory. In 
particular, Bacmann et al. (2000) utilized ISOCAM imag- 
ing data (~ 7/im), while Ward- Thompson, Andre, & Kirk 
(2002) utilized the ISOPHOT photometer at 90, 170, & 200 
pm. 

In Table 1, we reproduce the observed fluxes that we 
adopt in our study. We choose the 450, 850, and 1300 /im 
data from Shirley et al. (2000), as the results are from an 
array, rather than a single element bolometer. The 90, 170, 
& 200 /jm fluxes are taken from Ward- Thompson, Andre, 
& Kirk (2002). In both cases, the uncertainties in Table 1 
are the quoted statistical uncertainties added in quadrature 
with an assumed 30 per cent calibration uncertainty in the 
data. 

The spectral energy distribution (SED) can help con- 
strain the density distribution. However, the constraint is 
not unique (Butner et al. 1991; Men'shchikov & Henning 
1997; Doty & Palotti 2002). On the other hand, the spa- 
tial distribution of intensity on the sky in the form of maps 
provides significantly more data with which to constrain the 
models (Adams 1991; Ladd et al. 1991). As a result, we also 



Table 1. Fluxes toward L1544 



Wavelength 


Flux Density 


Uncertainty" 


ref 


(H 


(Jy) 






90 


< 0.57 6 


n/a 


1 


170 


14.9 6 


(0.2,4.5) = 4.5 


1 


200 


18.6 6 


(0.3,5.6) = 5.6 


1 


450 


17.4 C 


(6.7,3.5) = 7.5 


2 


850 


3.64 c 


(0.18,0.73) = 0.75 


2 


1300 


0.27 d 


(0.04,0.05) = 0.07 


2 



a (Statistical, 30% cal.) uncertainties added in quadrature 
b Flux density in a 150 arcsec aperature 
c Flux density in a 120 arcsec aperature 
d Flux density in a 40 arcsec aperature 

1 Ward-Thompson, Andre, & Kirk (2002) 

2 Shirley, Evans, Rawlings, & Grcgcrscn (2000) 

Table 2. Previously inferred conditions for L1544 



Parameter Value ref 

Optical depth at 1300^m (n 3 oo) 2 X 10~ a 3 

Optical depth at 200/im (T200) 0.06 5 

Density distribution BE™ or PL(m = 2) b 2 

ISRF(G^ Mp ) 0.6 6 

Bolometric Luminosity (Z/bol) [Lq] 1-0 ± 0.3 1 

Central Luminosity (X/ cen tral) [Lq] < 0-1 5 

Diameter (D) [x 1000 AU] ~ 17.8 4 

Central mass d (M120) [M ] O- 4 ^^ 1 

Total mass (M to t) [M Q ] 2.7 ± 0.7 2 

H 2 column (Nn 2 ) [10 22 cm" 2 ] 6-13 3,5 



a Bonner-Ebert sphere 

b Power-law of the form n(r) = no(ro/r) m 
c ISRF strength relative to standard MMP 
d In a 120 arcsec diameter beam 

1 Shirley, Evans, Rawlings, & Gregcrscn (2000) 

2 Evans, Rawlings, Shirley, & Mundy (2001) 

3 Ward-Thompson, Motte, & Andre (1999) 

4 Bacmann ct al. (2000) 

5 Ward-Thompson, Andre, & Kirk (2002) 

6 Young, Shirley, Evans, & Rawlings (2003) 

consider the 450 fj>m and 850 /im intensity maps produced 
by Shirley, et al. (2000). 

The data above were analyzed in detail in the presenting 
papers in order to constrain the structure and conditions 
associated with L1544. While these approaches were semi- 
analytic in nature and/or assumed spherical symmetry, they 
represent the norm in analyzing long wavelength emission 
from dust. Consequently, these results provide a useful first- 
order approximation to the underlying source conditions, 
and are presented in Table 2. 



3 MODEL 

We have constructed detailed, self-consistent, three- 
dimensional radiative transfer models through dust. The 
model utilizes a monte-carlo approach combined with an 



approximate lambda iteration to ensure true convergence 
even at high optical depths. The model has been tested 
against existing 1-D (Egan, Leung, & Spagna, 1988); and 
2-D (Spagna, Leung, & Egan 1991) codes with good suc- 
cess. 

Based upon the input parameters discussed below, we 
solve for the dust temperature and radiation field at each 
point in the model cloud with a typical resolution of ~ 10 15 
cm. The emergent radiation, after subtraction of the back- 
ground ISRF, is then convolved with the appropriate tele- 
scope beam (Shirley et al. 2000) for comparison with ob- 
servations. We do not simulate the chopping, as the chop 
direction and size varied with time during the observations. 
Evans et al. (2001) found that chopping affected the radial 
profile in the outer regions and were unable to determine an 
outer radius. While we will give an outer radius in this paper, 
it should be interpreted with caution. We have simulated 
the contours only out to radii that are not strongly affected 
by chopping, so conclusions about the shape of the density 
profile to that point should be safe. The approach of simulat- 
ing the observations has the advantage that we can directly 
compare models to observational output, while minimizing 
the need for imposing outside approximations/assumptions 
to the observational data. 

3.1 Input data 

We adopt a triaxial ellipsoidal density structure for our cloud 
models, of the form 



n(x,y,z) = n 




where n is the number density, no is the reference density, 
(x,y,z) defines the position, (a,b,c) are parameters speci- 
fying the shape of the ellipsoid, and m is the dust density 
distribution exponent discussed in the following subsection. 
This form is consistent with fact that the observed distribu- 
tion of projected axial ratios is well matched by randomly 
oriented, intrinsically prolate cores (Myers et al. 1991; Ry- 
den 1996). It is also consistent with the nonaxisymmetric 
evolution of magnetically subcritical cores (Nakamura & Li 
2002). We follow this work and the symmetry in the NW- 
SE direction, and reduce our density structure to a prolate 
spheroid having the long (x) axis in the NW-SE direction, 
and being symmetric in the y— and z— directions. In this 
case, we take b = c, and define the axis ratio to be q = a/b. 
Based upon the observed isophotes, we consider axis ratios 
in the range 1.2 < q < 5. 

We assume that the dust density follows a power-law, 
with exponent m. This significantly simplifies the source pa- 
rameterization while maintaining the ability to consider var- 
ious amounts of central condensation. The power-law param- 
eterization is commensurate with similar spherical models, 
as the spherical average density distribution for any ellipsoid 
given by equation (1) is n sp h C ricai oc r~ m , and is given by 
n sp hcricaiM = / ddn(x,y,z)/ J dd = n C(a,b,c,m)r~ m , 
where C(a, b, c, m) is a constant that depends only upon the 
geometry, and the density distribution. Furthermore, this 
approach is consistent with both observations (e.g. Evans et 
al. 2001) and theoretical (e.g. Nakamura & Li 2002; Curry 
& Stahler 2001; Ciolek & Basu 2000) models of evolution 
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of magnetized clouds. Based upon previous spherical mod- 
els, we consider density distribution exponents in the range 
1.0 < m < 3.0. 

The outer radius is the distance along the long (x—) 
axis at which the cloud structure terminates, and is half 
of the adopted long-axis diameter, D. At this point, the 
cloud is presumed to mesh with the surrounding, more dif- 
fuse, medium. We consider diameters in the range 0.06 < 
D(pc) < 0.2 (12000 < D(AU) < 40000). The inner radius 
defined as the distance from the center of the model space 
to the edge of the first non-central cell is r ln ~ 0.008 x D. 

The interstellar radiation field (ISRF) is the primary 
heat source for the grains in most of our models. We 
adopt the ISRF compiled by Mathis, Mezger, & Panagia 
(1983; hereafter MMP). The strength of the ISRF, Gmmp = 
(adopted ISRF) / (MMP ISRF), is in the range 0.3 < 
Gmmp < 3.0. We also include the potential effects of a lumi- 
nous central source by considering 10~ 6 < L*(Lq) < 10~ 2 . 

The absolute dust density (or mass) is constrained 
by the optical depth to source center along the long 
axis via the optical depth at 450 /jm, namely T450 = 
j n(x) (Q(450fj,m)iva 2 )dx. Here n(x) is the number density 
specified above, Q is the grain absorption efficiency at 450 
/jm discussed below, and a is the grain radius taken to be 
0.1 fim. Since the grains are small compared to the wave- 
length of the majority of the incident and re-emitted radi- 
ation (having A > l^im), the grain opacity per unit mass is 
independent of the grain size (Ossenkopf & Henning, 1994). 
Test models with different grain sizes confirm this indepen- 
dence of our results with adopted grain size. Finally, we 
consider optical depths in the range 0.01 < T45o,iong < 0.5. 

We adopt the dust opacities in column 5 of Table 1 
of Ossenkopf & Henning (1994). They have been successful 
in fitting observations of both high-mass (e.g., van der Tak 
et al. 1999, 2000; Mueller et al. 2002) and low-mass (e.g., 
Evans et al. 2001) star forming regions. These opacities are 
calculated for grains grown by coagulation and accretion of 
thin ice mantles for 10 5 years at a density of 10 6 cm -3 . These 
conditions should be applicable to the cold, dense cloud core 
of L1544. 



3.2 Grid of Models 

A guided search of parameter space yielded a preliminary 
best-fit (base) model as judged by the \ 2 deviation between 
the model and observations (see Section 4) . In order to con- 
firm the best fit, a relatively fine grid of over 400 models 
were constructed surrounding this best fit model to test the 
quality of fit and the geometry of the parameter space. From 
the base model, sub-grids were constructed by varying pairs 
of parameters: amount and distribution of dust (r, m), the 
strength of heat sources (Gmmp, L„), and the source ge- 
ometry (D, q). A best fit model from these sub-grids was 
identified as the new base model. The process was iterated 
multiple times to ensure the best fit was found. A final grid 
was then made surrounding the best fit model, with the 
comparison to observations discussed in Section 4. 
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4 COMPARISON WITH OBSERVATIONS 

In order to compare the model predictions with observa- 
tions, the simulated and observed data were fit using the 
reduced x 2 statistic. For the 450 /im and 850 fim sky maps, 
the major and minor axes of the isophote contours were 
compared, such that 



X450 — 



and 



2 

X850 



2 2 



=i j=i 



(•"450,i, j.obs 7*450, i, J, model ) 



450, i,j 



6 2 

EE 



('"850,i,j,obs — rs50,i,j, model) 



850, i,j 



(2) 



(3) 



Here i corresponds to the number of the contour level, j 
corresponds to the axis (major, minor), r is the radius of 
the i th contour along the j th axis, a is the uncertainty in 
the radius, and the leading fraction denotes the j^zt term 
with ./V450 = 4 and iVsso = 12 surface brightness contour 
data points (major and minor axes for each contour). The 
quality of the fit of the SED is measured similarly, via 



2 

XSED 



It, 



(-Ffe.obs — F k,: 



(4) 



Here k specifies the flux data point, is the flux of the k th 
datapoint, cr*, is the uncertainty in that flux, and the leading 
| denotes the j^rj term with Ased = 6 flux data points. 
The observed data are taken from Table 1. 

The overall quality of fit was then taken from the aver- 
age of the three reduced chi-squares (assuming that the 450 
/im, 850 /im, and SED results are independent) as 



2 

Xto 



= (X450 + X850 + XSEd)/3- 



(5) 



This choice for overall fit quality is advantaged by the fact 
that in general none of the individual \ 2, s dominate the oth- 
ers for good fits. Furthermore, the individual % 2 ' s are poten- 
tially more sensitive to different parameters. For example, 
the SED is dependent on the dust mass (e.g., Hildebrand 
1983; Doty & Leung 1994), while the surface brightness is 
more dependent upon the density and temperature distribu- 
tion (e.g. Shirley, Evans, & Rawlings 2002; Doty & Palotti 
2002). As a result, the fit to Xtot requires a simultaneously 
strong fit to all data. Unless otherwise noted, it is this value 
of chi-squared that is discussed below. For reference, our 
best fit model has Xtot = 0.5. 



4.1 Optical depth and density distribution 

The amount of dust and its distribution can be parameter- 
ized by r and m. In Fig. 1, we plot \ 2 contours for a grid of 
models in the T450 — m plane. The models run are denoted 
by the crosses. There exists a slight degeneracy between T450 
and m, in that a range of models yield \ 2 < 1- This is sim- 
ilar to the degeneracy in the r — m parameter space found 
by Doty & Palotti (2002) when fitting the SEDs of spherical 
models. However, in this case, the spatial distribution of in- 
tensity (i.e. the sky maps) allow the degeneracy to be more 
nearly broken. 

Based upon these results, we infer that the optical depth 
at 450 fj,m can be constrained to 0.06 < T45o,iong < 0.24. 




1.5 2.0 2.5 

Dust Distribution Exponent [ m ] 



3.0 



Figure 1. Contours of \ 2 fit for the models in the T450 — m plane. 
Note that while there is some degeneracy, the region is limited, 
enabling a constraint on both T450, and m. 



These values are, however, defined along the long-axis of the 
ellipsoid. The optical depth along the short axis is a factor of 
~ 4 lower, or in the range 0.015 < T450, short < 0.06. This cor- 
responds to optical depths of 2 x 10 -3 < T1300, short < 8xl0 -3 
at 1300^m, and 0.06 < r 20 o, short < 0.24 at 200pm. These 
ranges are consistent with those of previous observations 
and spherical modeling. This is to be somewhat expected, 
as at these wavelengths the source is transparent, so that 
nearly all grains can be seen. As a result, the total amount 
of energy absorbed from the ISRF plays a slightly larger 
role than the spatial symmetry, meaning that the spheri- 
cal model 'average optical depth' is an appropriate starting 
point. 

The results of Fig. 1 also suggest that the dust den- 
sity distribution can be well fit by a power law of the 
form in equation (1), with an exponent 1.75 < m < 2.25. 
This is generally consistent with the previous results from 
spherical modeling which suggest that a power law having 
m — 2 would fit the data (Evans et al. 2001), so long as a 
self-consistent temperature distribution is adopted. Again, 
this is somewhat expected, as the density powerlaw can be 
viewed as essentially fixing the relative amounts of 'warm' 
outer dust to 'cold' inner dust. Consequently, for a relatively 
regular source which is transparent at the wavelengths of 
interest, it appears that the spherical average density distri- 
bution provide a reasonable first-order estimate. Most strik- 
ingly, when a self-consistent, depth-dependent temperature 
distribution, and the source asphericity are taken into ac- 
count, the data are well- fit by a singular power-law density 
distribution. 

In order to probe the ability of a Bonner-Ebert-like den- 
sity power law to fit the observations, we have also consid- 
ered models with a flattened density distribution within a 
sphere of radius r cut . When we allow the density distribu- 
tion for r > r cut to vary (keeping m, and T450 constant), 
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Figure 2. Density profiles for L1544. The spherical average of the 
best-fit ellipsoidal distribution is given by the long-dashed lines. 
Also plotted are the corresponding best-fit Bonner-Ebert sphere 
(dotted line) and spherical power law (solid line) from Evans et 
al. (2001). 



we find that a central flattening for r < 0.002 pc = 400 
AU will produce adequate fits. On the other hand, when the 
density profile is unchanged for r > r cut (allowing T450 to 
vary), acceptable fits are found for r cut < 0.008 pc = 1600 
AU. This is roughly consistent with the results of Evans et 
al. (2001) in the sense that our density at 1600 AU is ~ 10 6 
cm -3 , the same as the central density for their best-fitting 
BE sphere. In order to place r cu t < 1600 AU in context, 
we note that the size of the ~ 15 arcsec beam at 850/^m is 
~ 0.01 pc = 2000 AU. As a result, while we cannot rule out 
a flattened Bonner-Ebert-like central density distribution, it 
is restricted to r < 1600 AU, and is essentially unresolved / 
at the resolution limit of the observations. 

The best-fit density distribution and spherical compara- 
tors are given in Fig. 2. The solid and dotted lines are the 
best-fit power law and Bonner-Ebert spherical density pro- 
files from Evans et al. (2001). The long-dashed line is the 
spherical average of our best-fit ellipsoidal profile (see Sect. 
3.1). Notice that, in general, our best-fit model represents a 
denser, smaller cloud than the previous spherical modeling. 
Furthermore, our best-fit model maintains an m = 2 profile, 
for r < -Rshort- However, once r > 7i s hort, the spherical av- 
erage density falls much faster with m ~ 4 — 5 due to the 
averaging of high density on-cloud and low-density off-cloud 
positions in taking the spherical average. Interestingly, this 
steep density fall-off is consistent with the results of Abergel 
et al. (1996) and Bacmann et al. (2000), even though along 
any given ray from the center of the cloud, m — 2. 

Finally, we note that our temperature distribution is 
roughly consistent with the self-consistent spherical model- 
ing of Evans et al. (2001), and the analytic work of Zucconi, 
Walmsley, & Galli (2001) for L1544. In particular, we find 
an unattenuated (outer) spherical average dust temperature 



2.Q r 1 1 T 1 * T 1 " 1 1 1 V 




10-6 1Q-5 1Q-4 1Q-3 

Luminosity [ L (solar Luminosity) ] 

Figure 3. Contours of fit \ 2 m the ISRF-L, plane. Note that 
there is no degeneracy between the parameters, and that signifi- 
cant limits can be set. 

of ~ 15 K, in keeping with the work of Zucconi et al. (2001), 
and an inner spherical average temperature of ~ 5.4 K. The 
inner temperature is consistent with the Zucconi results for 
a high central density (riH = 4.4 x 10 6 cm -3 , but approxi- 
mately 2 K lower than the Evans et al. (2001) and Zucconi et 
al. (2001) results for power-law density distribution adopted 
here. This has two causes. First, the central density for our 
power law model is closer to that of Zucconi et al. (2001). 
Also, while we our short axis is only slightly more opaque 
than the spherical models of Evans et al. (2001), our long 
axis is more opaque by a factor of ~ 2. The increased opacity 
leads to less radiation pentration, and thus lower tempera- 
tures in our models. 

4.2 ISRF and central luminosity 

The strength of the external and internal heat sources are 
parameterized by Gmmp and L*. In Fig. 3, we plot the \ 2 
contours in the Gmmp-£* plane. Again, the crosses signify 
the models considered. In this case, we can see that there is 
no degenercy between these two parameters (similar to re- 
sults for fitting SEDs, Lis et al. 2001). This is due to the fact 
that the size and spacing of the outer contours is mostly de- 
termined by the external heat source (the ISRF), while the 
inner contours are determined by the opacity at the wave- 
length of the dominant heating radiation and the luminosity 
of the central source. Consequently, each parameter probes 
a different region of the spatial intensity distribution maps. 

Based upon the results of Fig. 3, we can constrain the 
strength of the external heat source, namely the ISRF, to 
be 0.8 < Gmmp < 1.6. This is significant in that while it 
is consistent with the MMP ISRF, it is approximately a 
factor of 1.5-3 larger than the reduced field used by Evans 
et al. (2001), and the attenuated field adopted by Young 
et al. (2004) in their spherical modeling. This difference is 



6 S. D. Doty et al. 



due to the fact that for wavelengths which dominate the 
absorbed radiation (A < 25/im), the spherical clouds are 
effectively larger (more optical depth at larger radii). Thus, 
the stronger radiation field found here is consistent with the 
smaller surface area, mass, and higher optical depth in our 
model versus the equivalent spherical modeling. 

Last, it is useful to note that we have considered the 
Gmmp — D plane of fits, finding that Gmmp is relatively 
unaffected due to the small optical depth of the models. 
Together with the discussion above, this result points to the 
potential to infer local conditions when modeling combined 
spectral and spatial data. 

Of perhaps even greater interest, it is also possible to 
constrain the central luminosity. We find L, < 10 -3 Lq. 
This result is commensurate with the fact that no IRAS 
source exists at this position (Beichman et al. 1986), and 
with the assignation of this source as a preprotostellar ob- 
ject (Ward- Thompson, Scott, Hills, & Andre 1994). The 
dominant factor in constraining the luminosity is the sur- 
face brightness distribution. For the source mass and opti- 
cal depth, a central source of luminosity L > W~ 3 Lq yields 
a relatively strong point source in the maps, and a con- 
comitant decrease in the size and ellipticity of the intensity 
contours. 

The luminosity constraint at first appears inconsistent 
with the central luminosity quoted in Table 1. However, the 
limit in Table 1 arises from the difference between the inci- 
dent luminosity from the ISRF (~ O.ILq) and the emergent 
luminosity (~ 0.2Lq) found by Ward- Thompson, Andre, & 
Kirk (2002). As discussed by those authors, the small dif- 
ference suggests that L C entral ~ may well be within the 
uncertainties in their analysis. Consequently, the work here 
provides a much more strenuous constraint on the upper 
limit to the luminosity of any central source. Furthermore, 
while there may be infall (Tafalla et al. 1998; Williams et 
al. 1999), the low luminosity inferred here confirms the pre- 
vious result that the collapse is nearly isothermal with a 
correspondingly small compressional luminosity (Henriksen 
1994), and is inconsistent with the formation of a luminous 
protostar at the center of the collapsing core. 

4.3 Outer diameter and axis ratio 

The size and shape of the cloud are parameterized by D, 
and q. The \ 2 contours for the diameter (D) - axis ratio (q) 
plane are shown in Fig. 4. While there is some degeneracy 
between the diameter and axis ratio, it is again limited. In 
this case, the degeneracy is limited by the fact that the flux 
and hence the SED is more sensitive to the outer diameter, 
while the spatial intensity contours are sensitive to both the 
axis ratio and the diameter. 

Based upon the results in Fig. 4, the outer diame- 
ter of the long axis of the source can be constrained to 
0.08 < D(pc) < 0.16 (16000 < £>(AU) < 32000). This is 
consistent with the size inferred by Bacmann et al. (2000). It 
is however, smaller than the diameter adopted in the spher- 
ical modeling of Evans et al. (2001), who adopted an arbi- 
trarily large diameter to simulate chopping into an extended 
cloud. It should be noted that the actual cloud size may be 
somewhat larger than the value inferred here as discussed in 
Section 3. 

It is also possible to constrain the axis ratio of the 
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Figure 4. Contours of \ 2 fit f° r models in the D — q plane. Note 
that while there is a degeneracy, it is limited so that constraints 
can be placed on these parameters. 

source. We find 1.7 < q < 2.5. This is roughly consistent 
with the isophotes for the 450 fj,m and 850 /jm sky maps, 
which show axis ratios in the range of 1.5 — 4. The observed 
axis ratios are large because the external heating radiation 
must travel a longer physical distance along the short axis 
than along the long axis to reach the same optical depth. As 
a result, isophotes are more elongated than the underlying 
isodensity contours. 

4.4 Comparison of best-fit model with 
observations 

Based upon the results above, our best fit model for L1544 
is a prolate spheroid having an axis ratio of q = 2, with an 
ellipsoidal power-law density distribution given by equation 
(1) having m = 2. We infer the optical depth along the long- 
axis at 450 /im to be T450 = 0.12, and the central luminosity 
to be L, = 0. Finally, we infer that an unsealed MMP ISRF, 
and an outer diameter along the long axis of D = 0.1 pc = 
20000 AU (though it could be somewhat larger) provide a 
best fit to the data. 

A comparison of the 450 and 850 /im sky maps with 
the best fit model is shown in Fig. 6. Notice that the fits 
are generally quite good, in terms of the size and width of 
the intensity contours. For comparison, Fig. 7 gives similar 
sky-map intensity contours for representative models with 
Xtot — 2 (solid lines), and 5 (broken lines) respectively. The 
contour levels are the same as in Fig. 5. Notice that the 
sky maps for even these somewhat low values of Xtot are 
significantly worse than for the best fit model(s). 

Likewise, the observed and modeled SEDs are compared 
in Fig. 8. The observational data (symbols with error bars) 
are taken from Table 1. The solid line corresponds to the 
best-fit model, while the dotted and dashed lines correspond 
to the Xtot = 2, 5 models discussed above respectively. As 
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Figure 5. Intensity contours for L1544. Top panel: 450 fim; low- 
est contour level at 50% of peak flux (3 a), with levels increasing 
by 33%. Lower panel: 850 fim; lowest contour level at 20% of 
peak flux (3 a), with levels increasing by 13%. Data from Shirley, 
Evans, Rawlings, & Gregersen (2000). 



before, the best fit model reproduces the observed data well. 
On the other hand, although the Xtot = 5 model is well- 
outside the observed SED limits, the same is not true for 
the Xtot = 2 model. While the \tot = 2 model is at the very 
limits of acceptability, it is still consistent with the observed 
SED. This is in keeping with the result of Doty & Palotti 
(2002) that fitting the SED alone can yield degenerate solu- 
tions (see also, e.g., Men'shcikov & Henning 1997; Lis et al. 
2001). Furthermore, it underlines the significance of resolv- 
ing the source, and simultaneously fitting both the surface 
brightness and SED to best constrain the source properties. 
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Figure 6. Intensity contours for L1544 at for the best fit model. 
Top panel: 450 fim. Lower panel: 850 fim. The contour levels are 
the same as in Fig. 5. 



Finally, the inferred ranges and best fit values for the 
parameters specifying the structure and conditions of L1544 
are summarized in Table 3, along with values inferred from 
previous work. As noted previously, we see that conclusions 
from earlier modeling of observations at wavelengths for 
which the source is optically thin yield results which are 
relatively close to the best fit models. Even more interesting 
is the fact that under the same conditions, spherical mod- 
els yield results which are reasonable first-order approxima- 
tions. Still, it should be noted that even in the simple case 
of a centrally condensed, externally heated, optically thin, 
ellipsoidal power law density distribution with an axis ratio 
of ~ 2, geometrical and radiative transfer effects do cause 
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Figure 7. Model intensity contours for L1544 at 450/^m (top 
panel) and 850/im (bottom panel) for a representative model hav- 
ing Xtot = ^ (solid lines) and Xtot = ^ (dotted lines) . The contour 
levels are the same as those in Fig. 5. 



deviations from results inferred from simpler analyses. As 
a result, self-consistent, multi-dimensional modeling is an 
important tool for understanding radiation from aspherical 
sources. 



5 CONCLUSIONS 

We have constructed models for the non-spherical protostel- 
lar core L1544 utilizing a fully three-dimensional continuum 
radiative transfer model. After convolving the model output 
with the actual SCUBA beam profile, we are able to corn- 



Figure 8. Comparison of the SED for the best-fit model (solid 
line), a model with x 2 ot = 2 (dotted line), Xtot = ^ (dashed line), 
and observations (solid symbols with errorbars). 



Table 3. Inferred conditions for L1544 





Previous 


This 


Parameter 


Value 


Work 


Optical depth (T1300) [x10~ 3 ] 


2 


5 ± 3 


Optical depth (T200) 


0.06 


0.15 ±0.09 


Density power law (m) 


2 a 


2.0 ±0.25 


ISRF (Gmmp) 


0.6 


1.0 ±0.5 


Central luminosity (£*) [Lq] 


< 0.1 


< 10~ 3 


Diameter (D) [xlOOO AU] 


~ 17.8 


12 b - 24 c 


Central mass (Mf 20 ) [M ] 


4 +0 - 8 


0.2 ±0.1 


Total mass (M to t) [Mq] 


2.7 ± 0.7 


1.4 ±0.9 


H 2 column (JV H2 ) [10 22 cm" 2 ] 


6-13 


14 ±9 



a Bonner-Ebert sphere or Power-law with m = 2. 

b Diameter of short axis 

c Diameter of long axis 

d In a 120 arcsec diameter beam 



pare our results with existing observations, including both 
SEDs and sky maps. Based upon this work, we find that: 

1. It is possible to constrain the optical depth and den- 
sity distribution. While a degeneracy in this plane exists as 
noted by Doty & Palotti (2002), the spatial intensity distri- 
bution more nearly breaks the degeneracy. In particular, we 
find 2 x 10~ 3 < Ti3oo,short < 8 x 10~ 3 , and 0.06 < r 20 o,short < 
0.24. Likewise, we find that an ellipsoidal density power law 
with an exponent m ~ 2 (1.75 < m < 2.25), can fit the 
observations well. While a Bonner-Ebert-like centrally flat- 
tened density distribution cannot be ruled out, the central 
flattening can be restricted to r < 0.008 pc = 1600 AU 
(Sect. 4.1). 

2. The ISRF and central luminosity are nearly orthog- 
onal parameters. The ISRF can be constrained to be within 
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~ 50 per cent of that defined by Mathis, Mezger, & Pana- 
gia (1983), somewhat higher than previously inferred from 
spherical modeling. The central luminosity can also be con- 
strained to < 1CT 3 L Q , consistent with the assignation of 
this source as a pre-protostellar object (Sect. 4.2). 

3. The diameter and axis-ratio show a limited degener- 
acy, which is again broken by the spatial intensity distribu- 
tion. We infer a long-axis diameter of 0.08 < -D(pc) < 0.16 
(16000 < D(AU) < 32000), and an axis ratio 1.7 < q < 2.5. 
The outer diameter and axis ratio may each be somewhat 
larger due to potential on-source chopping in the observa- 
tions, and the projection of the long axis onto the plane of 
the sky. Both the diameter (D) and the axis ratio (q) are less 
than would be inferred directly from observations, and are 
due to radiative transfer effects and non-spherical geometry 
(Sect. 4.3). 

4. The values and ranges of values found here are 
roughly consistent with those inferred directly from observa- 
tions or spherical modeling, due to the relative transparency 
of the source at the wavelengths of observation. In general, 
we infer a smaller source of lower mass and higher optical 
depth / column density than previously assumed (Sect. 4.4). 

5. The SED yields degenerate fits. A "best-fit" requires 
both the SED and surface brightness distribution. This ap- 
proach, combined with a large grid of models (~ 400), can 
yield a best-fit and range of uncertainty in inferred parame- 
ters. Interestingly, a range of 10 (0.5 — 5) or even 4 (0.5 — 2) 
in x 2 can yield remarkably different quality of fit (Sect. 4.4). 



ACKNOWLEDGEMENTS 

We are grateful to Chad Young and Rebecca Metzler for 
interesting discussions. This work was partially supported 
under grants from The Research Corporation and the Deni- 
son University Research Foundation (SDD), NASA grants 
NAG5-10488 and NNG04GG24G to the University of Texas 
at Austin (NJE), and a Battelle Intership (MLP). 



Ladd, E. F., Adams, F. C, Fuller, G. A., et al. 1991, ApJ, 382, 
555L 

Lis, D. C, Serabyn, E., Zylka, R., & Li, Y. 2001, ApJ, 550, 761 
Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212 
Men'shcikov, A. B., & Henning, Th. 1997, A&A, 318, 879 
Mueller, K. E., Shirley, Y. L., Evans, N. J. II, & Jacobson, H. R. 

2002, ApJS, 143, 469 

Myers, P. C, Linke, R. A., & Benson, P. J. 1983, ApJ, 264, 517 
Myers, P. C, Fuller, G. A., Goodman, A. A., & Benson, P. J. 

1991, ApJ, 376, 561 
Nakamura, F., Li, & Z.-Y. 2002, ApJ, 566, 101L 
Ossenkopf, V., & Henning, Th. 1994, A&A, 291, 943 
Rydcn, B. S. 1996, ApJ, 471, 822 

Shirley, Y. L., Evans II, N. J., & Rawlings, J. M. C. 2002, ApJ, 
575, 337 

Shirley, Y. L., Evans II, N. J., Rawlings, J. C, & Gregersen, E. 

M. 2000, ApJS, 131, 249 
Spagna, G. F. Jr., Leung, C. M., & Egan, M. P. 1991, ApJ, 379, 

232 

Tafalla, M., Mardones, D., Myers, P. C, et al. 1998, ApJ, 504, 
900 

van dcr Tak, F. F. S., et al. 1999, ApJ, 522, 991 
van dcr Tak, F. F. S., et al. 2000, ApJ, 537, 283 
Ward-Thompson, D., Andre, P., & Kirk, J. M. 2002, MNRAS, 
329, 257 

Ward-Thompson, D., Motte, F., & Andre, P. 1999, MNRAS, 305, 
143 

Ward-Thompson, D., Scott, P. F., Hills, R. E., & Andre, P. 1994, 

MNRAS, 268, 276 
Williams, J. P., Myers, P. C, Wilner, D. J., & di Francesco, J. 

1999, ApJ, 513, L61 
Young, C. H., Shirley, Y. L., Evans, N. J., & Rawlings, J. M. C. 

2003, ApJS, 145, 111 

Young, K. E., Lee, J.-E., Evans II, N. J., Goldsmith, P. F., & 

Doty, S. D. 2004, ApJ, in press 
Zucconi, A., Walmsley, C. M., & Galli, D. 2001, A&A, 650, 662 



REFERENCES 

Abergel, A. et al. 1996, A&A, 315, L329 
Adams, F. C. 1991, ApJ, 382, 544 

Andre, P., Ward-Thompson, D., & Motte, F. 1996, A&A, 314, 
625 

Bacmann, A., et al. 2000, A&A, 361, 555 

Beichman, C. A., et al. 1986, ApJ, 307, 337 

Benson, P. J., & Myers, P. C. 1989, ApJS, 71, 89 

Butner, H. M., Evans II, N. J., Lester, D. F., Lcvreault, R. M., 

& Strom, S. E. 1991, ApJ, 376, 636 
Ciolek, G. E., & Basu, S. 2000, ApJ, 529, 925 
Curry, C. L., & Stahler, S. W. 2001, ApJ, 555, 160 
Doty, S. D., & Leung, C. M. 1994, ApJ, 424, 729 
Doty, S. D., & Palotti, M. L. 2002, MNRAS, 335, 993 
Egan, M. P., Leung, C. M., & Spagna, G. F., Jr. 1988, Comput. 

Phys. Comm., 48, 271 
Elias, J. H. 1978, ApJ, 224, 857 

Evans II, N. J., Rawlings, J. M. C, Shirley, Y. L., & Mundy, L. 

G. 2001, ApJ, 557, 193 
Henriksen, R. 1994, in Montmerle T., Lada C. J., Mirabel I. F., 

Tran Than Van J., eds., The Cold Universe, Editions Fron- 

tieres, Gif-sur-Yvette, p. 241 
Hildebrand, R. H. 1983, QJRAS, 24, 267 



